AD-A213  494 


NEW  TRANSIENT  ALGORITHMS  FOR 
NON-NEWTONIAN  FLOWS 


David  S.  Malkus, 
Yi-Cheng  Tsai, 
and 

Robert  W.  Kolkka 


UNIVERSITY 
OF  WISCONSIN 


CENTER  FOR  THE 
MATHEMATICAL 
SCIENCES 


Center  for  the  Mathematical  Sciences 
University  of  Wisconsin— Madison 
610  Walnut  Street 
Madison,  Wisconsin  53705 


October  1989 


(Received  October  11,  1989; 


Sponsored  by  Approved  for  public  release 

Distribution  unflmHed 

U .  s.  Arm/  Research  Office 
p.  0.  Box  12211 
Research  Triangle  Park 
North  Carolina  27709 


Air  Force  Or  rice  or 

Scientific  Research 
Washington,  DC  20332 


National  Science  Foundation 
Washington ,  DC  20550 


89  TO  ?0  2  54 


UNIVERSITY  OF  WISCONSIN-MADISON 
CENTER  FOR  THE  MATHEMATICAL  SCIENCES 

NEW  TRANSIENT  ALGORITHMS  FOR  NON-NEWTONIAN  FLOWS  * 

David  S.  Malkus1 
Yi-Cheng  Tsai1 
Robert  W.  Kolkka  2 

CMS  Technical  Summary  Report  #90-14 
October  1989 

Abstract 

A  fully  dynamic  method  for  shear  flows  is  presented  that  treats  the  short  time-scales 
associated  with  Newtonian  viscosity  (or  short  relaxation  processes)  and  shear-wave  prop¬ 
agation  implicitly,  while  treating  the  long  relaxation  processes  explicitly.  The  method  is 
generalized  to  flows  with  non-constant  strain-rate  histories  in  the  context  of  the  well-known 
fiber-drawing  problem.  The  linearized  stability  of  the  methods  is  analyzed,  and  extension 
of  these  methods  to  planar  flows  is  given.  The  approach  taken  in  the  case  of  non-trivial 
deformation  histories  is  that  of  an  Oldroyd  difference  quotient  (ODQ)  that  approximates 
the  convected  derivatives  of  the  differential  constitutive  equation  in  Lagrangian  fashion 
along  the  portion  of  the  streamline  upstream  of  the  stress  evaluation  point.  Techniques 
based  on  earlier  ideas  of  drift-function  tracking  are  used  to  develop  a  weighting  scheme  for 
the  ODQ  that  permits  the  use  of  low  order,  Co  stress  elements.  The  numerical  methods 
are  discussed  and  analyzed  in  the  context  of  a  Johnson- Segalman  fluid  model  with  added 
Newtonian  viscosity.  The  resulting  initial-boundary-value  problem  is  globally  well-posed 
and  possesses  the  key  feature:  the  steady  shear  stress  is  a  non-monotone  function  of  the 
strain  rate.  Such  models  will  be  seen  to  display  the  spurt  phenomenon  in  plane  Poiseuille 
flow  and  apparently  related  phenomena  in  step  strain  experiments.-'  Analysis  of  the  nu¬ 
merical  methods  shows  that  the  ratio  of  short  to  long  relaxation  times  or  the  Newtonian 
viscosity  ratio  is  a  key  parameter  in  the  stability  and  accuracy  of  the  methods.  When  this 
is  properly  accounted  for  the  techniques  described  here  work  well  in  shear  and  extensional 
flows  and  show  promise  for  two-dimensional  flows.  "  > 
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1.  INTRODUCTION 


We  consider  the  equations  of  motion,  stress,  and  continuity  for  a  Johnson-Segalman 
fluid  [1]  with  added  Newtonian  viscosity.  The  Johnson-Segalman  model  should  be  viewed 
as  a  simple  constitutive  equation  representing  features  which  are  characteristic  of  non¬ 
monotone  stress/strain-rate  constitutive  relations  [2  -  10]  (see  Fig.  1);  many  other  models 
have  these  features,  and  some  are  known  to  produce  the  kind  of  dramatic  and  interesting 
rheological  behavior  discussed  in  this  paper  in  the  context  of  Johnson-Segalman  [2,4,6]. 
There  is  some  controversy  among  rheologists  as  to  whether  the  predictions  of  such  models 
are  physically  correct;  we  believe  so,  and  have  argued  our  case  extensively  elsewhere  [2.4  - 
7,10].  Our  thrust  here  is  numerical,  and  we  aim  to  develop  methods  of  general  applicability 
to  models  with  and  without  non-monotone  behavior.  At  very  least,  the  success  of  our 
methods  in  shearing  flows  in  the  ‘spurt’  regime  induced  by  non-monotonicity  can  be  viewed 
as  pushing  the  numerical  method  to  the  limit  in  order  to  test  its  robustness  in  the  presence 
of  strong  instabilities. 

The  equations  of  motion,  stress,  and  continuity  in  terms  of  velocity,  v,  total  stress,  S, 
extra  stress,  r,  strain  rate,  e,  and  pressure,  p,  for  the  Johnson-Segalman  model  are 

[&+(v-V)v]=V-S 
S  =  p  I  -f-  2se  +  r 

<  %f+r  =  2^  (1) 

%f  :=  It  +(V'V)r-rVv  -  (Vvf  r  +  (1  -  a)  ( re  +  er) 
k  V  •  v  =  0  . 

The  equations  are  in  nondimcnsional  form;  a  is  a  ratio  of  Reynolds  number  to  Deborah 
number:  time  has  been  scaled  by  the  dominant  relaxation  time;  e  is  the  ratio  of  Newto¬ 
nian  viscosity  to  polymer  zero  shear  viscosity  or,  alternatively,  the  ratio  of  dominant  to 
secondary  relaxation  time:  the  dimensional  form  of  the  equations  and  further  details  of 
tin'  dimensional  analysis  may  be  found  in  Refs.  [2]  and  [5].  The  parameter,  a,  is  a  nondi- 
mensional  number  that  relates  the  motion  of  the  polymer  molecules  to  the  mean  motion 
of  the  continuum:  this  motion  is  non- affine  when  a  /  1:  we  consider  only  a  >  0  here.  In 
Refs.  [5.6.7],  it  is  shown  that  the  Newtonian  viscosity  term  is  generally  representative  of 
the  presence  of  shorter  relaxation  times,  widely  separated  from  the  dominant  relaxation 
time  that  can  occur  in  pure  polymer  systems  as  well  as  polymer/solvent  systems. 

In  Refs.  [5.6.7],  the  reduction  of  Eqs.  (1)  to  transient  shear  flows  is  discussed.  Fol¬ 
lowing  Ref.  [5],  for  a  <  1  one  may  also  scale  a.  v,  and  /  by  (1  —  a2)1/2;  the  result  is 

'  or,  =  <7r  4-  £v i z  4-  / 

<  at  =  (Z  4-  1)  vt  —  a  (2) 

k  Z,  =  -avt  -  Z  . 

In  the  absence  of  Newtonian  viscosity  (i.e.,  when  z  =  0),  the  system  (2)  is  hyperbolic  as 
long  as  Z  4-  1  >  0.  and  the  speeds,  c,  are 


c  =  0,  ± 


Z  4-  1 


Of 


(3) 


9 


lo9,oTi2 


lo910  "Y 

Fig.  1  Stress  vs.  strain-rate  for  a  Johnson-Segalman  model.  Di¬ 
mensional  units:  values  characteristic  of  actual  materials.  Values 
of  s  give  neutrally  stable  (1/8)  and  critical  cases  (0.015);  7  =  2ej2. 

Otherwise'  the'  system  has  elliptic  character,  and  the  type  change  corresponds  to  a  loss  of 
rvoluf ionarity.  When  Newtonian  viscosity  is  present  (i.e.,  when  s  ^  0),  it  can  be  shown 
that  the  system  is  always  evolutionary  [S.9].  When  Newtonian  viscosity  is  absent  the 
system  can  exhibit  Hadainard  instabilities  [S.9]. 

2.  NEW  TRANSIENT  ALGORITHMS  FOR  SHEAR  FLOW 
Tn  this  section  we  briefly  describe  two  test  problems  that  are  solved  by  our  numerical 
techniques  for  shear  flow.  We  then  describe  the  numerical  method,  arrange  the  algorithm 
in  suitable  form  for  a  linearized  'matrix'  stability  analysis  [11.12],  and  derive  the  resulting 
stability  bounds. 

,4.  Spurt:  In  Refs.  [2,4  -  7],  the  phenomenon  of  spurt  observed  by  Vinogradov,  et  al.  [13] 
is  modeled  exploiting  the  non-monotonicity  of  the  Johnson-Segalman  models  steady  stress 
vs.  strain-rate  curve,  using  system  (2)  and  boundary  conditions  appropriate  for  pressure- 
driven  How  (sec  Fig.  ?).  The  phenomenon  of  the  dramatic  increase  in  capillary  or  die 
throughput  at  a  fixed  critical  stress,  independent  of  molecular  weight,  can  be  reproduced 
with  remarkable  accuracy.  The  resulting  flow  regime  is  characterized  by  a  layer  of  high 
shear  rate  near  the  wall,  in  which  the  extra  stress  plays  virtually  no  role  in  equilibrating  the 
driving  pressure  gradient.  Outside  the  layer,  the  longer  relaxation  process  is  the  dominant 
load-bearing  mechanism.  During  the  process  of  spurt,  there  are  large  oscillations  in  stress 
and  shear  rate,  and  material  points  that  end  up  inside  the  layer  are  alternately  subjected 
to  high  and  low  rates  of  shear  until  effective  steadiness  is  achieved. 

D.  Step  strain.*:  We  consider  an  idealized  model  of  the  classical  step  strain  experiment 
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x  =  -h/2 


x  =  0 

Fig.  2.  Problem  domain  for  channel  How  driven  by  a  pressure 
gradient  (plane  Poiseuille  How).  The  dimension,  h/2,  is  scaled  to 
1/2  in  the  non  dimensional  equations.  T  =  a-\-evx  is  the  total  shear 
stress. 


[3,10,14,15]  that  combines  Eqs.  (2)  with  boundary  conditions 


(  v(0,t)  =  0 

<u(M)  =  0(<)  (4) 

{  <rx(0,t)  -f-  evxx(0,t)  =  0  . 

We  thus  have  Couette  flow  [3j  in  which  the  upper  wall  is  a  plate  that  moves  suddenly, 
shortly  after  time  t  =  0.  The  picture  is  similar  to  Fig.  2,  with  the  symmetry  line  replaced 
by  the  fixed  wall,  whose  coordinate  is  taken  to  be  x  =  0  in  nondimensional  units;  the 
moving  wall  is  at  x  =  1  (see  Fig.  5,  below).  The  classical  problem  is  properly  posed  [10] 
by  taking  the  Unlit  of  the  sequence  of  solutions  to  this  problem  as  g  tends  to  a  multiple 
of  the  Heaviside  function,  which  is  the  step  in  strain.  7,  based  on  plate  motion.  In  this 
context,  we  can  see  the  conseqviences  of  adding  an  arbitrarily  small  amount  of  Newtonian 
viscosity,  c,  or  omitting  this  term.  We  consider  the  Lodge-Meissner  function: 


$(7)  =  lim 

f  — *00 


Til  ~  r22 

7  S12 


2 £ 

lT  • 


(5) 


where  T  :=  a  +  srx,  and  7  is  the  strain  calculated  from  system  (3),  which  is  a  factor 
of  (1  —  a2)1//2  smaller  than  the  corresponding  strain  of  system  (1)  (the  a  =  1  case  is 
treated  as  a  special  case).  Note  that  the  total  shear  stress,  Si 2  or  T,  appears  in  the 
denominator  of  Eq.  (5).  Eq.  (5)  was  intended  to  apply  in  situations  where  there  is  no  flow 
with  relaxing  stresses  and  thus  where  there  is  no  difference  between  S\ 2  and  772  (or  T  and 
<7):  here  we  have  chosen  to  use  the  physically  measurable,  total  stress,  for  reasons  that 
will  become  clear  later.  The  Lodge-Meissner  relation  [10,14,15]  states  that  the  limit  of  <E> 
as  the  step  becomes  instantaneous  is  observed  to  be  1  in  real  polymer  systems.  For  the 
Johnson- Segalman  model,  Eq.  (5)  can  be  evaluated  analytically,  under  the  assumption  that 
momentum  effects  can  be  ignored  (even  though  the  step  is  assumed  to  be  instantaneous) 
and  the  flow  remains  homogeneous  for  ail  time  [10].  The  result  is  plotted  in  Fig.  3.  The 
Jolmson-Segalman  model  has  been  criticized  for  the  obvious  failure  to  satisfy  the  Lodge- 
IvieiWit-r  relation  and  the  fact  that,  the  shear  stress  evidently  changes  sign  as  a  function  of 
strain,  as  evidenced  by  the  pole  in  Fig.  3. 


C.  The  numerical  method  for  shear  flow :  In  these  shear-flow  examples,  Newtonian  viscos¬ 
ity  dominates  in  high  shear-rate  regions  or  times;  there  the  Eqs.  (2)  have  the  parabolic 
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Fig.  3.  The  Lodge-Meissner  function  with  a  =  0.95.  Values  for 
7  <  5  are  close  to  1.  The  pole  is  at  7  «  10.  Another  pole  is  at 
7  ~  30;  poles  periodically  repeat. 

character  of  Newtonian  flow;  implicit  treatment  of  Newtonian  viscosity  is  essential.  In  the 
regions  or  times  of  lower  shear-rate,  where  the  longer  relaxation  time  dominates,  there  is 
wave  propagation  at  the  elastic  wave  speed  given  by  Eq.  (3),  and  the  equations,  though 
not  formally  classifiable,  behave  like  a  hyperbolic  system  with  a  small  amount  of  added 
viscous  damping.  For  high  polymers  of  the  type  used  by  Vinogradov,  et  al.  [13],  the  wave 
speed  is  finite  but  very  high,  and  numerical  methods  must  respect  it  by  either  the  time- 
step  restrictions  of  explicit  integration  or  the  cost  of  implicit  integration.  On  the  other 
hand,  phenomena  such  as  the  latent  development  of  spatial  inhomogeneity  in  step  strains 
take  place  in  several  step  rise-times;  phenomena  like  Vinogradov’s  spurt  take  time  on  the 
order  of  many  dominant  relaxation  times  to  unfold.  In  either  case,  the  time  scale  of  the 
interesting  physical  phenomenon  is  orders  of  magnitude  longer  than  wave  traversal  times. 
To  capture  such  phenomena  by  explicit  methods,  controlled  by  element  wave  traversal 
times,  is  totally  out  of  the  question.  The  following  method  was  developed  to  confront  this 
difficulty  without  requiring  implicit  iteration  on  the  nonlinearities  [2,5]. 

—  Spatial  discretization:  We  use  standard  finite  element  interpolations  [11,12].  In  what 
follows,  the  definitions  of  shape  functions,  matrices,  etc.  are  at  the  global  (rather  than 
element)  level.  The  shape  functions  are  linear,  one  basis  function  per  node,  as  picture  in 
Fig.  4.  Following  standard  practice,  the  nodal  values,  V{,  are  functions  of  time  alone,  and 
the  shape  functions.  A7,,  are  functions  of  space  alone: 


[N1(x),N2(x),...,Nxi(x)}  =  [N] 

{vi(t)  \ 

I  =  (A'l  {*>}  • 

VM{t)  ) 


(6) 
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Fig.  4.  Velocity  trial/ test  functions  for  Galerkin  method.  The  Nj 
arc  basis  functions  for  piecewise  linear  approximations  for  velocity. 
Stresses  arc  assumed  to  be  piecewise  constant  on  cells  between 
nodes;  cells  have  length  L}  and  may  van ■  in  size. 


We  use  standard  "B-Matrix"  notation  for  the  matrix  relating  strain  rates  to  nodal  values. 
Thus 


Vr  =  [Z?]  {l’}  =  [Nl,t,N2,r . iV.U, :]{v}  . 

The  acceleration  can  be  written  as  follows: 


(7) 


=  [-V] 


r  ) 


=  [.V]  {v}  . 


(S) 


'  l'M.t  f 

The  viscoelastic  contribution  to  momentum  diffusion  is 


f‘  fl 

/  <pTadx  =  /  {6}T  [B]Tcrd.c  , 

J  0  Jo 


(0) 


in  which  vve  shall  assume  a  is  elementwise  constant.  The  Galerkin  form  of  the  momentum 
equation  is 
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a<f>vtdx  +  /  4>xadx+£  /  <pzvzdx  =  0  , 


r** 

Jo 


for  all  admissible  d,  at  each  t  >  0.  When  the  admissible  test  functions,  <f>,  are  expanded 
in  terms  of  the  A',,  the  spatial  discretization  of  momentum  equation  becomes  [11,12] 


[M]{a}+e[A']M  +  {PF}  =  {0}  , 


where 


M  =  =  {«} 


[M}=  f  a[A:]T[Ar]  dx 
Jo 

{W}  =  f  [D]r a  dx 

Jo 

[/,']=  f[B)T[B)dx 

Jo 


[A']  =  [A'1(.r),iV2(x),...,WA/(x)] 

[5]  =  [Ari.*,A2.r . ArA/,x]  • 

—  Time  differencing:  The  time  differencing  is  based  on  a  standard  generalized  trapezoidal 
scheme  [11.12]. 

{  rn  +  l  }  —  {l’n}  +  At  [(1  —  C){an}  +  C{an+1  }]  ,  (13) 

where  £  >  0  is  a  parameter  of  the  method.  A  predictor  for  the  viscoelastic  contribution  to 
tin-  shear  stress  (the  extra  stress)  is  given  by  a  mixed  difference,  forward  on  the  stress  and 
nonlinear  terms  and  backward  on  the  linear  strain-rate  term.  Z  and  a  are  assumed  to  be 
constant  on  each  element  and  advanced  bv 


° n +1  & n 


—  Zr i(ez)n  +(t’x)J1+[  dn  . 


The  momentum  equation  is  solved  based  on  this  prediction,  and  a  corrector/Z-updato 
cycle  is  used  for  both  stress  equations. 


d  n-fl  d  n 

At 

Z  n  +  1  ~  Zn 

At 


—  ( Zn  +  1 )  C’x )  n+  i  ® n 
=  —  d „  + 1  ( vz  )ri^_  j  -  Z n  . 


Combining  the  space  and  time  discretization  of  the  momentum  equation,  we  hav< 

[A/]{«n  +  .}+£[A’]{W„  +  ,  }  +  {»*„+,  }  =  {0}. 


(1") 


where  {W„+i}  is  a  predictor  of  the  viscoelastic  contribution, 

{hFn+1}  =  f  [B]Tan+i  dx  .  (18) 

Jo 

With  the  discretization  of  system  (2),  whose  stability  is  analyzed  below,  corrector  is 
not  used,  and  there  are  no  “~”ed  corrected  quantities.  We  thus  have 

- -  =  Zn(vx)n  +  (ux)n+l  ~  an  ,  (19) 

and  Z  given  by  Eq.  (16)  (with  dn+i  replaced  by  crn+1).  It  is  possible  to  apply  the  a- 
corrector  in  the  frozen  Z  case,  but  numerical  experiments  show  that  this  procedure  is 
less  stable  than  the  one  analyzed  here.  Furthermore,  experiments  also  show  the  nonlinear 
algorithm  with  the  predictor-corrector  scheme  “nearly”  unconditionally  stable  when  pa¬ 
rameters  are  selected  according  to  frozen  Z  analysis,  without  the  corrector  step.  We  should 
einphcisize  that  we  are  analyzing  a  somewhat  simpler  algorithm  -  without  the  correction 
cycle  -  to  approximate  the  behavior  of  the  more  complicated  algorithm.  Experience  shows 
that  this  is  a  valid  simplification. 

—  Computational  arrangement(momentum):  When  we  use  Eq.  (19)  to  predict  crn+\  in 
Eq.  (17),  the  result  is 


{TFn+i}=  [  [B]r<rn+ldx 
Jo 

=  At\K){vn+])  +  At  j  Zn[B]T{B}{vn}  dx  +  (l  -  At)  (  [B]Tandx 

Jo  Jo 


Substitution  into  momentum  equation  yields 


[d/]{u„-(-t }  +  (£  +  Zlt)[A]{t’n+i}  +  {IVn+1 }  —  {0}  , 


vhere 


Z  —  term 


{I Vn  +  l}:=At  /  Zn[B}T[B}{vn}dx  +  (l-  At) 


f  [B]T a n  i 

Jo 


Note  that  the  for  frozen  coefficient  system  (2),  the  Z-term  in  {Wn+i}  can  be  written  as 
At.Z[K}{vn). 

—  A  recursion  for  {U'n+1 };  From  Eq.  (20)  and  the  definition  of  {IFn+i},  we  have 


/W  <7  n+1  dx  =  At[K}{vn+l)  +  {Wn- 
Jo 


Cycling  back  to  n+1  — +  n  yields 


f  [B]TCr n  dx  =  Zto[A']{Vn}  +  {W'n}  • 
Jo 
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We  substitute  into  definition  of  {irn+i}  to  get 


{H'n+i}  =  (1  -  4i){H*n}  4-  4t(l  -  At)[K}{vn}  +  At  f  Zn[B)T[B){vn)  dx  .  (25) 

J  0 

In  the  frozen  Z  case,  Eq.  (25)  becomes 

{tfn+i}  =  (1  -  +  (Z  +  (l~  At))At[I\}{vn}  .  (26) 

—  Quasi  linear  difference  equations:  The  fully  nonlinear  system  (1)  can  now  be  written  in 
recursive  form. 

MK+i}  +  (e  +  At)[K}{vn+l}  +  {HVh}  =  {0} 

{ti'n+i}  =  (l-At){Wn)+At(l-At)[K]{vn}  +  At  [  Zn[B]r[B]{e„}dx  (27) 

Jo 

{t'n  +  l}  -  C^{«n  +  l}  =  {Vn}  +(1  -  C)-U{«n}  , 

after  each  cycle  of  which.  Zn  +  l  is  obtained  from  Eq.  (16)  (<rn+1  =  crn+1);  though  it  is  not 
necessary  for  algorithmic  closure,  an+1  can  be  obtained  from  Eq.  (19)  after  each  cycle. 
The  frozen  coefficient  form  of  Eq.  (27)  is 

MK+i }  +  (=  +  &t)[K]{vn+l }  +  {U'n+1 }  =  {0} 

{irn+!}  =  (1  -  At){\vn}  +  At(l  -  At)[K]{vn}  +  AtZ[K}{v„}  (2S) 


{  *’n+  1  }  —  C-lt{«n  +  l  }  =  {t’n}  +  (1  —  (  )~lt  {nn  }  . 

Note  that,  for  the  current  scheme  <r„ ,  Zn  and  {vt)n  =  [B]{u„}  are  all  piecewise  constant 
in  .r. 

—  Modal  decomposition:  Let  ['!>]  be  a  matrix  of  generalized  eigenvectors  of  [A']  -  sf.U], 
normalized  so  that 


WhA'IM  =  [!)]  =  ding(a.-,2 ) 

(^riA/im  =  [i\ . 

We  will  investigate  the  stability  of  the  frozen  coefficient  linearized  system.  Let 

{«n  +  i}  =  ['L]{dn  +  1} 

{!’n+l}  =  [T]{0n  +  1}  . 

Tiie  momentum  equation  gives 


(29) 


(30) 
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where 


{«»+i }  +  (^  +  ^i)[fi]{rn+i }  +  {II  (i4-i }  =  {0} 


(31) 


{U'-„+1}  =  (»I-]7'{U;+1}  .  (32) 

Suppressing  the  subscript  V’  and  writing  a  typical  component  of  {dn+1}  as  simply  a„+1. 
etc.  The  decoupled  system  is 

^n  +  I  +  (i  +  +  l  n+  1  =0 

II  =  (1  —  At) IVn  +  At[  1  —  At)uj^  vn  +  AtZuP'v-n  (33) 

t'ri+1  C*-^^h+1  —  t-?n  T  (1  C)Atan  . 

Though  =  0  is  ruled  out  by  boundary  conditions  in  the  problems  considered  here,  the 
lowest  frequencies  in  the  mesh  are  physically  meaningful,  and  the  u  -*  0  limit  must  be 
considered.  It  is  easy  to  see  with  our  simple  elements  [11,12]  that  the  highest  frequencies 
are  mesh-dependent,  and 


ui 


2 

max 


max 


nip 


(34) 


D.  Stability  <uuili/.<i.<:  d he  goals  of  stability  amdysis  are  twofold:  first,  to  give  a  prior 
estimate  of  the  critical  time  step,  if  any.  and  second,  to  assure  that  the  numerical  method 
does  not  have  the  effect  of  artificially  suppressing  or  delaying  the  onset  of  the  instability. 
The  onset  of  spurt  corresponds  to  the  instability  of  the  frozen-coefficient  linearization.  We 
will  then  analyze  the  discrete  system  (33)  applied  to  the  same  frozen-coefficient  system 
and  find  that  the  instability  in  tin*  discrete  system  occurs  at  precisely  the  same  values  of 
the  parameters  as  the  exact  linearized  system.  The  exact  solution  can  be  expanded  in 
Fourier  series:  a  typical  mode  is 


c(.r,/)  =  f(t)  shikar  <r(.r.  t )  =  ij(t)  cos*c.r  .  (35) 

Substituting  into  system  (2) 


«/'(0  =  --~,2/(0  -  <]{fU 

g'(t)  =  (Z  +  \Uf(t)-g(t)  . 


The  system  eigenvalues  are 


£a»’*/n  -f-  1  / ( tu,'*/o  4-  1 )“  (Z  - f-  1  4-t)u,'^ 

- 2 - ±  V - 5 - o - 


(37) 


be  s,  ri  t  <  0  only  if  Z  +  1  -ft  >  0  (decaying  solutions)  —  the  Equation  Stability  Condition. 
The  phase-plane  analysis  of  Refs.  [G.7]  shows  that  this  condition  is  the  same  condition  that 
characterizes  the  end  of  the  'latency  period  and  the  onset  of  the  spurt  instability.  We 
shall  derive  stability  bounds  assuming  that  u,’  £  (O.+oc). 
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—  The  stress  relaxation  condition:  The  recursion  for  U'n  is  based  on  the  underlying  O.D.E. 
for  cr,  which  in  time-discrete  form  is 

an+l  —  (1  —  At)(7n  +  At(Zn(vx)n  +  (l’T)n+i  )  ,  (3S) 

which,  when  there  is  no  flow,  is 


aI1+1  =  (1  -  At)an  .  (39) 

This  is  a  discretization  of  (exponential)  stress  relaxation,  the  most  fundamental  non- 
Xewtonian  phenomenon.  It  is  no  restriction  to  assume  that  At  is  always  small  enough 
for  accurate’  stress  relaxation.  The  non-oscillatory  stress  relaxation  condition  (X’OSRC)  is 
tln’ii. 


At  <  1  .  (40 j 

The  XOSRC'  is  necessary  to  accurately  reproduce  the  lowest  frequency  phenomena. 
To  analyze  the  stability  of  higher  frequencies,  we  write  the  first  equation  of  Eqs.  (33)  in 
step  n  and  n  -f  1;  solving  this  for  nn,  an 4.,  and  substituting  into  the  remaining  equations 
of  the  system  (33).  results  in  tli  following  2x2  system: 


where 


6  AW2(Z  +  6)  1 


A]  - 

nc-it+d >;ArA{Z+9)+J-\ 

L  A  A 

B  =  0(At)  =  (1  -  At) 

1  - 

A  = 

1  -  ( AtiS 

i  =  £  +  At  . 

We  can  bound  the  spectral  radius  by  using  the  invariants  of  [A]  (see  [11.12] ): 


(41) 


1 42) 


"i  =  otr[Al  =  0 

1 


[  X  +  B)  +  3  —  1 

"  + - 5 - 


(43  1 


a,  =  det [A]  =  —  [d(J-  1)  -  (i  -  C )At2J2(Z  +  B 
Spi’ctral  radius  of  [A]  being  <  1  is  characterized  by  [11.12] 

Cl  2  4“  1  <7  ■>  1 

-  — -  <  a!  <  — - ,  a  >  <  1  , 

—  1  —  rv  '  * 


’44' 


and  algebraic  growth  is  ruled  out  by 


—  1  <  a  1  <  1.  aj  =  1  . 


(45) 
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—  Extreme  frequencies:  It  turns  out  that  by  considering  the  u  — >  0  and  u  — *  oo  behavior 
of  Eqs.  (43)  -  (45),  the  most  salient  features  of  the  method  described  here  can  be  identified; 
a  more  detailed  analysis  is  given  elsewhere  [16].  When  w-»0we  find  that 


The  roots  of  equation, 


«.-=(«  +  !). 


a  2 


6  . 


o 

S  —  2diS  +  02  =  0  , 


(46) 


are  B ,  and  1;  thus  requiring  \B\  <  1  implies  At  <  2.  This  condition  is  irrelevant,  since  the 
NOSRC  supercedes. 

When  u  — »  oc, 


P  Too, 


P~i  ,  (1-0 

A  C 


and 


{Z  +  B)At2u2  ( Z  +  B)At 

A  * 

SO 


(Z  +  B)At  1-C' 

£  C  J  ’ 

6(  1-C)  ,  (l-c )(Z  +  d)At 

a  >  -> - 7 - H - t - r -  • 

C  C  £ 

Thus  the  quadratic  Eq.  (46)  has  roots 

(47) 

For  (  <  1/2  conditional  stability  is  anticipated,  and  specific  analysis  for  large  u>  is 
required.  This  analysis  is  carried  out  in  detail  in  Ref.  [16],  where  it  is  found  that  there 
is  a  dependence  of  At  on  Ax2  that  is  characteristic  of  parabolic  problems  [11,12],  arising 
from  the  Newtonian  viscosity  term.  Since  the  algorithm  is  semi-implicit  and  requires  the 
factorization  of  matrices,  there  is  no  gain  in  taking  C  <  1/2.  In  Ref.  [16],  it  is  also  shown 
that  for  C  >  1/2.  there  are  no  interior  values  of  the  spectral  radius  between  u.’  =  0  and 
u.*  =  oc  higher  than  those  bounded  by  the  NOSRC  and  by  requiring  that  the  roots  of 
Eq.  (47)  be  less  than  or  equal  to  1.  This  implies  the  following  stability  summary  for 
(>1/2: 

( 1 )  From  —  0  4-  ( Z  4-  9)At/i  <  1,  we  deduce  that  (t  4 -  Z  —  1)_R  —  2i  <  0.  If  £  +  Z  —  1  <  0. 
there  is  no  condition.  Otherwise 


(1-C)  ni(Z  +  B)At 

s  =  - ; — ,  —u  -j - r - 

C  £ 


At  < 


2s 

s  4-  Z  —  1 


12 


(2)  From  —d  +  (Z  +  6)At/z  >  —1,  we  deduce  that  (Z  +  1  +  z)A t  >  0.  If  Z  +  1  -f  e  >  0. 

there  is  no  condition.  Otherwise  the  numerical  method  and  the  linearized  equations 

are  unstable. 

From  the  phase  plane  analysis  of  Refs.  [5,6,7],  it  can  be  deduced  that  for  the  range 
of  parameter  space  appropriate  for  modeling  Vinogradov’s  polyisoprenes  [13],  Z  is  always 
negative  unless  it  is  forced  to  be  positive  by  artificial  initial  conditions.  Furthermore,  Z  <  0 
whenever  the  primary  normal  stress  difference,  jVj,  is  positive,  as  it  is  for  any  steady  flow 
of  a  Johnson-Segalman  fluid.  Circumstances  in  which  Case  (1)  above  would  hold  require 
that  Ni  be  substantially  negative  and  have  not  been  observed  by  the  authors.  However, 
the  linearized  equations  do  not  account  for  the  Z-equation,  thus  could  apply  equally  well 
to  a  shear  thickening  fluid  in  which  Z  >  0  is  likely;  in  such  cases,  the  condition  of  Case  (1) 
could  apply.  In  Case  (2),  whenever  the  Equation  Stability  Condition  holds,  the  algorithm 
is  essentially  unconditionally  stable  for  all  time  steps  smaller  than  that  implied  by  the 
NOSRC,  which  is  an  accuracy  condition,  independent  of  Ax.  When  the  exact  equations 
are  unstable,  the  numerical  method  also  reproduces  solutions  that  grow  rather  than  decay 
in  time.  Thus  the  onset  of  the  spurt  instability  is  not  suppressed  by  an  artificial,  numerical 
increase  in  the  value  of  Z  +  1  +  z. 

The  key  to  the  success  of  the  algorithm  lies  in  the  first  equation  of  system  (27)  that 
predicts  the  extra  stress  using  the  zero  shear  viscosity  implicitly  and  the  non-Newtonian 
contribution,  represented  by  the  nonlinear  Z  term,  explicitly.  No  nonlinear  equations  are 
solved  when  time  is  advanced  according  to  the  second  equation  of  system  (27)  ,  but  the 
fastest,  wave  speed  is  that  at  zero  shear,  as  can  be  seen  from  Eq.  (3)  and  the  fact  that 
Z  <  0.  Then'  is  no  restriction  on  the  time  step  in  the  linear  stability  analysis  of  Eqs.  (6)  - 
(13),  and  in  practice,  time  steps  as  large  as  a  half  million  element  wave  traversals  can  be 
stably  employed.  However,  the  analysis  of  Refs.  [5,6,7]  shows  that,  when  z  is  small,  the 
shear-flow  problem  is  extremely  stiff:  there  is  a  characteristic  response  time  of  the  system 
of  0(z )  (the  ‘‘Newtonian  phase”)  in  which  the  polymer  contribution  to  the  shear  stress,  rr. 
develops.  During  this  short  period,  Z  changes  substantially  but  then  continues  to  change 
on  a  time  scale  that  is  typically  much  longer,  determined  by  the  dominant  relaxation  time, 
which  is  scaled  to  1  in  the  nondimensional  equations.  In  order  for  the  slow  change  in  Z  to 
be  accurately  resolved,  the  short-time  change  in  Z  must  have  been  accurately  resolved.  In 
Ref.  [5]  it  is  shown  that  the  stiffness  requires  time  steps  on  the  order  of  fractions  of  £  at 
early  times  to  obtain  accurate  solutions.  Time  steps  can  be  lengthened  somewhat  during  a 
quiescent,  period  ( "latency" ),  but.  as  shown  in  Refs.  [6.7],  the  0(s)  scale  reemerges  in  the 
dynamic  spurt  process.  The  point  is  that  the  semi-implicitness  can  be  used  to  accurately 
smooth  out  wave  propagation  at  early  time,  since  elastic  waves  play  no  important  role  in 
the  spurt  phenomenon;  however,  this  leaves  an  intricate  interplay  between  two  additional, 
widely  separated,  time  scales  that  must  be  properly  accounted  for  in  numerical  solutions. 

3.  NUMERICAL  RESULTS  FOR  STEP  STRAINS  IN  SHEAR  FLOW 

To  summarize  what  occurs  in  numerical  simulations:  when  z  ^  0,  the  total  shear 
stress,  T.  has  a  near  ^-function  contribution,  and  for  t  >  0  there  is  spatially  homogeneous 
stress  relaxation  when  a  =  1.  Spatial  inhomogeneity  develops  for  a  <  1  (at  sufficiently 
large  7  [10]).  With  e  =  0,  the  stresses  are  bounded  for  .all  time  if  a  =  1,  but  there  is 
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FIXED  WALL 


t=0 . 50  t=0 .  60  t=0. 70 

KOVING  WALL 


FIXED  WALL 


Fig.  5.  Velocity  profile  as  time  evolves  for  dynamic  solution  of 
system  (2)  with  boundary  conditions  (4).  a  =  O(10~12).  a  =  0.95. 
c  =  0.001417.  and  a  step  of  duration  0.1  see.  At  this  average  strain 
of  7  =  12.  the  solution  develops  a  spatial  inhomogeneity  in  a  latent 
spurt  with  flow  in  the  opposite  direction  from  plate  motion.  The 
fixed  wait  i*  at  j  —  0  and  the  moving  wall  at  x  =  1. 


Hadamard  instability  if  a  <  1  and  7  is  sufficiently  large.  Thus  the  effect  of  Newtonian 
viscosity  is  significant.  For  a  =  1,  its  presence  does  not  make  a  difference  in  the  satisfaction 
of  the  Lodge-Meissner  relation,  but  it  does  make  an  increasingly  large  difference  in  the  total 
force'  required  to  move  the  plate  as  the  steps  approach  instantaneous  rise.  For  a  <  1.  the 
Lodge-Meissner  relation  is  not  satisfied  at  large  strain,  but  in  the  presence  of  Newtonian 
viscosity,  the  computed  solutions  at  large  strain  arc  not  necessarily  unphysical,  as  has 
been  claimed  [15].  As  is  described  in  Ref.  [10],  the  failure  to  satisfy  the  relation  can  bo 
associated  with  an  inhomogeneous  flow  regime  that  develops  some  time  after  plate  motion 
has  ceased  (Figs.  5-7). 

Apparently  a  similar  phenomenon  to  spurt  occurs  in  flows  of  a  Johnson-Segalman 
fluid  after  a  sudden  step  in  strain.  This  is  illustrated  in  Fig.  5:  There  is  an  initial  period 
in  which  an  essentially  homogeneous  flow  occurs.  For  smaller  strains,  roughly  7  <  5  for 
data  corresponding  to  Vinogradov,  et,  aids  polyisoprenes,  the  homogeneity  peisists  for 
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Fig.  6.  Local  strain  at  the  walls  vs.  time  for  dynamic  solution  of 
system  (2),  with  boundary  conditions  (4),  a  =  0(1O~12),  a  —  0.95, 
i  =  0.001417,  and  a  step  of  duration  0.1  sec.  At  this  average  strain 
of  7  =  S,  the  solution  is  inhomogeneous  at  long  times. 


all  time.  For  moderate  strains,  roughly  5  <  7  <  10  for  Vinogradov's  data,  a  spurt -like 
inhomogeneous  flow  develops  while  the  stresses  are  relaxing  and  the  flow  has  essentially 
ceased.  At  moderate  strains,  the  spurt  is  in  the  direction  of  plate  motion.  For  large 
strains,  roughly  7  >  10.  the  spurt  is  in  the  opposite  direction  from  the  plate  motion,  as  it 
is  in  Fig.  5.  These  observations  explain  the  inhomogeneity  in  local  strain  based  on  fluid 
motion  of  Figs.  G  and  7  (in  contrast  to  7,  an  average  strain,  based  on  plate  motion).  It  is 
interesting  to  note  that,  the  formula  of  Fig.  3  is  derived  by  ignoring  momentum  effects  that 
are  responsible  for  the  development  of  spatial  inhomogeneity,  but  the  position  of  the  pole 
in  that  formula  correctly  predicts  the  direction  of  the  spurt.  These  phenomena  are  not 
fully  understood  at  the  time  of  this  writing,  but  it  seems  clear  that  the  criticisms  of  the 
Johnson-Segalman  model  based  on  assumptions  of  spatial  homogeneity  are  too  simplistic. 

4.  HIGHER  DIMENSIONAL  FLOWS  -  FIBER  DRAWING  AND  2-D 

To  generalize  these  ideas  to  multidimensional  flow  requires  approximation  of  the  stress 
gradients  in  the  converted  derivative  of  Eqs.  (1).  These  can  be  avoided  by  observing  that 
the  effective  strain  tensor  [3( p.  4SG )]  of  the  Johnson-Segalman  model  satisfies 

f  E,(/')  =  — E,(f')(ne-u,)(0  (4S) 

\  E((f)  =  I  , 
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strain 


Fig.  7.  Local  strain  at  the  walls  vs.  time  for  dynamic  solution  of 
system  (2)  with  boundary'  conditions  (4),  a  =  O(10-12),  a  =  0.95, 
e  =  0.001417,  and  a  step  of  duration  0.1  sec.  At  this  average  strain 
of  7  =  12,  the  solution  is  inhomogeneous  at  long  times,  but  the 
inhomogeneity  is  reversed  from  what  it  is  for  strains  on  the  other 
side  of  the  pole  in  the  Lodge-Meissner  function  (see  Fig.  3). 


where  t1  represents  historical  time,  and  u  is  the  spin  tensor.  A  Lagrangian  difference 
quotient  along  particle  pathlincs  can  be  used  to  approximate  the  whole  convected  derivative 
of  system  (1)  [3(p.  494)].  and  the  time-discrete  stress  equations  become 


where 


Var  r"+1  -  Et+At(t)TnEj+^t(t) 


Vt 


At 


=  2en+1  -  rn 


(49) 


rn  :=  r(x(f).f) ,  r"  :=  r (x(f  +  At),t) ,  rn+1  :=  r (x(f  +  At ), t  4-  At ))  .  (50) 

Since  this  formulation  represents  an  approximation  to  the  Oldroyd  derivative  [17]  as  a 
difference  quotient,  we  refer  to  it  as  an  “Oldroyd  Difference  Quotient”  (ODQ  hereafter). 
Thus  the  formulation  of  the  simple  shear  problem  can  be  employed  with  piecewise  constant 
stresses,  since  the  full  stress  gradient  field  is  not  required.  Existing  particle  tracking 
methods  can  be  used  to  evaluate  the  required  strain  tensors  [IS]  with  the  velocity  field 
held  fixed  between  time  steps. 
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The  key  to  the  success  of  the  ODQ  formulation  of  Eqs.  (50)  is  the  calculation  of  an 
appropriate  value  of  rn.  Eqs.  (49)  reduce  the  evaluation  of  the  convected  derivative  to 
the  one-dimensional  calculation  of  an  ordinary  derivative  along  a  streamline,  but  along 
that  streamline,  the  basic  stress  approximations  are  piecewise  constant.  However,  the 
problem  of  extracting  derivative  information  from  a  piecewise  constant  approximation  to  a 
differentiable  function  of  one  variable  seems  preferable  to  the  employment  of  2D,  Co  stress 
fields  and  the  evaluation  of  (weak)  stress  gradients.  Here  we  propose  scheme  based  on 
a  transit  time-weighted  average  between  the  value  of  the  stress  in  the  element  containing 
x(t  +  At)  and  the  value  in  the  adjacent  element  upwind  of  x(t  + At)  on  the  same  streamline. 
The  weighting  is  tuned  to  recover  O(At)  accuracy  in  Eqs.  (49).  Approaches  to  2D  flows 
that  incorporate  these  features  immediately  suggest  themselves;  here  we  illustrate  our 
approach  in  the  fiber  drawing  problem  of  Fig.  9  [19,20].  The  mathematical  and  physical 
ramifications  of  this  simplified  model  of  fiber  drawing  have  been  thoroughly  investigated 
elsewhere  (see  Refs.  [19  -  21],  and  additional  references  cited  in  Ref.  [22]  for  details); 
furthermore,  the  finely-tuned  methods  of  Refs.  [19]  and  [20]  seem  to  be  the  methods  of 
choice  for  this  problem.  Our  purpose  in  solving  this  problem  is  to  illustrate  the  potential 
of  the  ODQ  formulation  in  flows  with  non-trivial  history  dependence  (i.e.,  the  substantial 
derivative  does  not  reduce  to  an  ordinary  one)  in  a  model  problem  which  can  be  posed 
in  one  spatial  dimension,  and  which  results  in  a  relatively  small  system  of  PDEs.  The 
problem  will  also  serve  to  reemphasize  the  importance  of  proper  numerical  treatment  of 
multiple  time  scales  in  a  rather  different  context  from  the  shear-flow  problem. 


extruder 


swell 


x,v 


Fig.  8.  The  domain  and  coordinate  system  for  the  fiber  drawing 
problem. 


A.  Mathematical  formulation:  We  employ  the  same  nondimensionalization  as  used  for 
Eqs.  (1)  and  (2);  the  nondimensional  domain  is  the  interval  (0,1),  and  time  is  scaled  by 
the  relaxation  time  [2,5],  etc.  In  the  problem  domain,  the  polymer  contribution  to  the 
stress  tensor  (the  extra  stress)  is  assumed  to  be  of  the  form  [21] 
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Fig.  0.  Staggered  grid  for  fiber  drawing  problem.  Velocity,  v  and 
area,  .4,  arc  at  grid  points,  while  stresses,  r  and  and  velocity- 
gradient,  vz,  are  at  half-points.  Problem  makes  no  sense  if  velocity- 
changes  sign  in  threadline;  positive  v  assumed. 


*  r  0  O' 

r  =  0  "  0  (51 ) 

_0  O  u 

The  velocity  field  is  irrotational  with  gradients  given  by 

0  0 

K  0  .  (52) 

0  -\vz. 

The  total  drawing  (x-direction)  stress  consists  of  extra  stress,  pressure  (radial  stress),  and 
Newtonian  contributions: 


Vv  =  2e  =  0  - 

0 


T  :=  N  +  3svz 


(53) 


where 


N  :=  t  —  zr  . 


(54) 


In  terms  of  these  quantities  and  the  area  of  the  fiber,  A{x,t),  the  equations  of  motion,  con¬ 
tinuity  and  stress  reduce  to  a  4x4  system  [21],  which,  combined  with  boundary  conditions 
is 


IS 


'  oc A(vt  +  vvx )  =  (TA)X 
At  -f-  ( vA)x  =  0 

*  Tt  +  VTX  -  2avx  =  -T  +  2vx  (55) 

+  vzox  +  avx  -  —  w  —  vx 

,  u(0,f)  =  De,  4(0,1)  =  1,  r(0,<)  =  ro,  vx(0,t)  =  vxO. 

Our  concern  here  is  solution  procedures  of  the  stress  equations  and  not  the  momentum 
equation,  so  it  is  particularly  appropriate  to  eliminate  the  latter,  making  the  same  reduc¬ 
tion  appropriate  to  Vinogradov’s  data  [13],  which  leads  to  the  system  of  ODEs  analyzed  in 
the  phase  plane  in  Refs.  [5,6,7].  When  a  is  formally  set  to  zero  in  fiber  drawing,  however, 
the  result  is  a  3x3  system  of  PDEs 


< 


(  At  =  -vAx  -  vzA 
rt  -  -vtx  -  (1  -  2 avx)r  +  2vx 
zct  =  —vzax  —  (1  +  avx)za  —  vx 

f  /  A  —  T  +  ~ 


0  =  -Ur  + 


3c- 


(56) 


The  fourth  equation  is  an  algebraic  relation  that  can  be  substituted  into  the  first  three 
equations  and  be  integrated  with  respect  to  x  to  give  a  closed,  3x3,  hyperbolic  system  with 
integral  coefficients.  The  unknowns  arc  two  stresses  and  the  area;  the  three  wave  speeds 
an'  r. 


D.  Finite  difference  discretization:  We  envisage  that  future  planar  and  axisymmetric  algo- 
rit hms  implementing  the  ODQ  will  employ  finite  element  discretization  of  the  momentum 
equation:  in  the  current  reduced  problem,  it  is  straight-forward  and  instructive  to  use  fi¬ 
nite  differences  for  all  equations  of  system  (55).  Motivated  by  the  finite  element  shear-flow 
formulation,  wo  use  the  staggered  grid  scheme  pictured  in  Fig.  9.  This  makes  stresses  and 
strain  rates  cell  quantities,  while  area  is  a  grid  quantity.  The  algebraic  relation  between 
strain  rate  and  velocity  gradient  is  approximated  by 


(r*)n 


l  •— 


//.4"  -  r"  ,  +  , 


3c 


U"  :=  \(Al  +A]'n_l)  , 

while  a  midpoint  rule  is  used  to  obtain  grid  values  of  velocity. 


t-r  =  dc 

u"  =  v”  4-  A xj^  M?+ 1/2  m  >  1 


m  —  1 


i=l 


(5S) 


An  explicit  forward  in  time/backward  in  space  difference  scheme  is  used  for  the  area 
equation. 
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I  At 

(«*)"  :=  I 


An  -  An 


=  —  vr 


m  — 1 


Ax 


-(V:)nAnm 


(u*)”  1  +  ( vxT  l 


m+2 


2  J 


(59) 


Turning  to  the  ODQ  formulation  for  the  stress  equations,  we  begin  by  observing  that 
the  transit  time  between  two  points,  x  and  y,  on  the  threadline  in  a  steady  velocity  field, 
v(x,t)  =  u(x)  is  given  by  [23] 


6{X-V)  =  !SW)-  (60) 

This  is  a  simple  example  of  a  “drift  function”  [18];  it  will  be  used  to  define  the  appropriate 
components  of  fn.  Since  we  intend  to  apply  the  ODQ  to  finite  element  schemes  in  the 
future,  it  is  appropriate  to  use  finite  element  interpolations  of  the  kind  defined  using  [N] 
in  Sec.  2.C  and  pictured  in  Fig.  4  to  provide  interpolate  velocities  between  grid  points.  We 
further  approximate  velocities  as  being  constant  between  time  steps,  so  that  the  steady 
velocity  field  u(x)  in  Eq.  (60)  will  be  completely  specified  by  nodal/grid  values.  Referring 
to  Fig.  9,  we  find  that 


Infi  Ax 

Si  :=  6(td,xm- 1)  =  — ^ : ^ - 

vm  -  Vm  _  i 


(61) 


and 


62  <5(xm- 1 1 £u)  — 


i"0(1  +  ^ef))Ax 


^m- 2  vm  —  1 

We  approximate  the  logarithms  in  Eq.  (62)  by  a  Taylor  series  and  find  that 


(62) 


6 £ 


(63) 


where  i  and  j  are  m,  m  -  1,  or  m  -  2  as  indicated  by  Eq.  (62).  Note  that  =  62  = 
A.r/2vm-i  the  transit  time  from  £u  to  is  approximated  by  6C  :=  61  +  62  =  Ax/vm-i\ 
this  will  be  used  to  define  the  appropriate  components  of  fn . 

We  now  turn  our  attention  to  approximation  of  the  strain.  Let 


A  t 

k  :=  J  vz  (x(t'),  t')  dt'  . 
W<’  use  a  first-order  approximation  to  the  integral: 


A tvx(Zd,t) 

Si vx(£d,t)  +  {At  -  6,K«„,0 


6X  >  At 
Si  <  At  . 


(64) 


(65) 
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The  two  cases  represent  the  situation  in  which  the  particle  at  the  centroid  of  a  given 
element  stays  in  the  same  element  during  one  time  step  and  the  case  in  which  it  passes 
into  the  next  element  upstream,  respectively. 

Note  that  k  =  O(At)  during  a  single  time  step,  and  thus  a  Taylor  series  expansion 
for  the  strains  of  Eq.  (4S)  is  appropriate.  The  appropriate  strains  and  their  Taylor  series 
approximations  are  then 


„  f  -E,2,  =  exp(2a/c)  =  1  +  2 qk  s  =  t 

E  F  \  ~  i  W 

I  i^22  =  exp(  —  QK)  =  1  —  QK  s  =  xc  . 


The  following  definitions  will  allow  us  to  specify  the  ODQ  algorithm  for  a  generic  stress 
equation: 


r 


Ax 


2  s  =  r 

—  1  s  =  zu 


(67) 


Note  that  \b\  =  0(1)  with  respect  to  At,  and  b  =  0  in  the  corotational  [3],  or  a  =  0  case. 
In  terms  of  6,  a  generic  strain,  E  is 


E2  =  1  +  bAt  .  (68) 

We  define  a  generic  stress  (either  r  or  zu),  ■s^_1y2,  at  time-level  n  and  staggered  grid  point 
m  -  1/2.  The  appropriate  generic  component  of  fn,  sn,  is  defined  using  r,  which  can  be 
seen  to  be  the  fraction  of  the  transit  time  from  one  element  centroid  to  the  next  that  is 
covered  in  one  time  step.  Stability  restrictions  discussed  in  the  next  section  essentially 
assure  that  r  <  1.  s"  is  defined  by 

:=  (1  _  r)^m  — 1  /2  +  rsm- 3/2  *  (^) 

W’e  use  Eq.  (69)  and  take  the  two  stress  equations  in  the  form  of  Eq.  (55);  the  ODQ 
method  for  a  generic  stress  equation  is 


sn+\  -  (1  +  bAt) 

(1  -  r)sn  !  +  rsn  3 

m-2 

m  2  m  2  . 

At 


—  —Sn  1  -I-  c(vx)n  !  . 


(70) 


Note  that  the  whole  left  side  of  the  stress  equations  in  Eq.  (55)  is  approximated  by  the 
ODQ.  For  the  purposes  of  analysis,  we  rearrange  Eq.  (70)  as  follows: 


,n+l 


At 


ZU,  -  _„n  . 

~  vm- 1 


Ax 


—  S 


-v^_1bAt- 


Ax 


-(1  -b)8”  X+C(vz)n  , 


added  term 


(71) 
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The  ODQ  method  in  this  simple  problem,  with  the  approximations  we  have  made  for 
transit  times  and  strains  is  a  modified  backward  difference  scheme  along  the  streamline. 
It  reduces  to  a  simple  backward  difference  scheme  in  the  corotational  case. 


C.  Stability  analysis:  For  the  purposes  of  stability  analysis,  we  first  study  the  corotational 
case  where  a  =  0.  It  is  not  clear  whether  this  case  is  physically  relevant  to  actual  polymer 
systems,  but  the  equations  are  more  tractable,  and  the  results  of  the  analysis  of  this  case 
admit  to  a  straight-forward  modification  for  a  ^  0  that  works  well  in  practice.  In  the  a  —  0 
case,  the  discretization  of  the  area  equation  and  ODQ  treatment  of  the  stress  equations 
reduces  to  a  backward  difference  method  for  the  whole  system.  Since  the  wave  speeds  are 
all  v,  inspection  of  Eqs.  (56)  suggests  that  a  stability  bound  of  At  <  Ax/v  is  sufficient. 
This  is  because,  as  will  be  illustrated  in  the  subsequent  analysis,  the  principal  part  of  the 
system  usually  dominates  the  stability  bound,  while  the  source  terms  contribute  the  effect 
of  an  change  of  O(Ax)  in  v  in  the  conventional  bound.  Since  the  conventional  bound 
is  slightly  pessimistic,  a  stable  and  accurate  time  step  can  usually  be  obtained  by  using 
a  modest  safety  factor.  This  is  not  the  case  here  when  £  is  small.  First,  consider  the 
linearized  area  equation, 


S  4 

<5.4,  =  —v6Az  +  ~N  + 
06 


a  r 

3i  Jr, 


At  +  w)  dx'  +  Ajy  , 


(72) 


where  the  varied  quantities  are  preceded  by  “<5,”  and  the  expansion  point  by  corresponding 
quantities  without  8 s.  The  four  terms  on  the  right-hand  side  have  the  following  description: 
principal  part  (highest  order  term),  intermediate  order  source  term,  low-order  (integral) 
source  term,  and  coupling  term  (intermediate  order),  respectively.  First  we  consider  the 
behavior  of  the  linearized  area  equation  uncoupled  from  the  stress  equations  by  taking 
8X  —  0.  A  standard  Neumann  analysis  [24,25]  shows  that  the  backward  difference  method 
applied  to  the  linearized  area  equation  is  stable  in  the  sense  required  for  convergence 
[25(p.  42)]  if  the?  conventional  bound  on  At  is  enforced.  This  restriction,  however,  allows 
unbounded  growth  for  t  — ►  co  on  a  fixed  grid,  due  to  growth  of  low  frequency  disturbances. 
This  arises  from  the  fact  that  the  intermediate  order  term  (with  coefficient  proportional 
to  .V)  is  typically  positive.  Numerical  experiments  with  the  backward  difference  method 
applied  to  Eq.  (72)  show  that  the  integral  term  is  not  sufficient  to  eliminate  the  growth 
of  low  frequency  disturbances.  Such  unbounded  growth  can  be  ruled  out  by  the  coupling 
terms,  if  the  linearized  stress  equations  governing  8 X  are  stably  integrated.  We  turn  our 
attention  to  a  criterion,  based  on  the  stress  equations  alone,  that  will  prevent  blow  up 
of  discrete  solutions  as  t  — ►  oc.  We  pursue  this  course  since  we  are  interested  in  finding 
steady  solutions  on  fixed  grids. 

A  simple  bound  can  be  derived  by  further  uncoupling  the  stress  equations  fr  n  each 
other.  This  bound  turns  out  not  to  be  sufficient,  but  it  illustrates  the  essence  of  the 
problem.  We  linearize  a  typical  stress  equation  and  set  the  perturbed  quantities  from 
other  equations  to  zero,  and,  as  in  the  case  of  the  area  equation,  we  ignore  the  variation 
of  the  velocity  that  is  reflected  in  the  lowest  order  (integral)  term.  Letting  q  :=  and 
considering  a  disturbance  of  wave  number  uj/k  in  a  standard  Neumann  analysis,  yields  a 
growth  function,  p ,  given  by 
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where  v  is  the  velocity  field  of  the  expansion  point  of  the  linearization.  We  note  that  |pj 
is  bounded  by  1  for  to  =  0  if  At  <  2/  ^1  -f-  ;  this  restriction  is  always  superseded  by 

the  bound  arising  from  cjAx  =  ±tt  (wave  number  1/Ax),  which  turns  out  to  be  the  most 
unstable  mode.  This  is 


At  < 


(74) 


The  bound  is  most  stringent  where  v  is  maximum,  usually  at  x  =  1. 

There  seems  to  be  no  accuracy  requirement  to  make  Ax  <C  e\  good  results  will  be 
exhibited  below  in  cases  where  quite  the  opposite  is  true.  Yet  in  such  cases,  the  bound  of 
Eq.  (74)  is  much  smaller  than  the  conventional  bound,  and  this  cannot  be  avoided  in  this 
explicit  formulation.  What  we  see  is  another  manifestation  of  the  extreme  stiffness  induced 
by  a  small  amount  of  Newtonian  viscosity;  in  the  semi-implicit  shear-flow  algorithm,  it 
only  affects  accuracy  during  those  periods  of  dynamic  response  when  the  “Newtonian" 
time  scale  is  active.  In  the  explicit  fiber-drawing  algorithm,  it  affects  stability,  and  the 
Newtonian  time  scale  must  be  respected  at  all  times.  Since  this  stiffness  arises  from  the 
competition  of  the  wave  speeds  of  the  principal  part  with  the  time  scale  inherent  in  the 
source  terms,  we  have  dubbed  it  “source  term  stiffness." 

Numerical  experiments  show  that  the  bound  of  Eq.  (74)  is  over-optimistic  by  roughly 
a  factor  of  3/2.  We  write  the  coupled  stress  equations  in  matrix  form. 


{t-;?,/,}  =  {«•;.■/,}  -  <•«  ({c-,/*}  -  {r:-,*})  -  §  w  {c,;_1/2} . 


(75) 


when’ 


[*]  = -3*- m  +  [  \ 


( 7G ) 


The  discrete  system  (75)  is  diagonalizable.  and  the  result  is  an  uncoupled  system: 
standard  Neumann  analysis  can  be  carried  out.  giving  a  bound  analogous  to  Eq.  (74).  in 
which  the  1  4-  -jj  in  the  denominator  of  Eq.  (74)  is  replaced  by  the  magnitude  of  the  most 
negative  eigenvalue  of  T[£?j.  The  result  is 


A  t  < 


Ax 


='  +  (!  +  ;)t  ' 


(T7> 


We  note  that  this  is  the  same  bound  that  would  result  subtracting  the  two  stress 
equations  to  produce  a  single’  equation  for  ,Y.  The  second  eigenvalue  of  [8]  is  0(1)  (with 
respect  to  c)  and  is  irrelevant  to  stability.  The  same  trick  of  combining  the  equations 
would  not  work  in  the  case  a  ^  0,  so  that  the  matrix  analysis  is  necessary.  We  linearize’ 
Eq.  (71):  as  can  be  seen  from  the  bound  of  Eq.  (77),  and  as  will  be  evident  in  the  numerical 
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results  presented  below,  the  crucial  terms  for  stability  analysis  are  those  of  order  0(s~ 1  ). 
The  only  terms  in  the  discrete  equations  (71)  with  a  ^  0  that  differ  from  our  previous 
analysis  that  we  chose  to  include  are  chose  of  order  0(c-1).  These  include  the  terms  of 
order  0(£-1)  already  accounted  for  in  the  analysis  of  the  corotational  case  and  additional 
terms  arising  from  “6b."  Furthermore,  among  those  terms,  those  that  arise  from  linearizing 
the  “added  term’’  in  Eq.  (71)  are  0(Af)  smaller  than  those  obtained  from  linearizing  the 
source  term  (1  —  so  we  ignore  the  former  and  include  only  the  latter.  We  note 

that  this  admittedly  ad-hoc  procedure  is  equivalent  to  retaining  only  0(e-1)  terms  in  the 
linearization  of  Eqs.  (56)  and  applying  a  backward  difference  scheme  (without  the  “added 
term")  to  the  result;  this  has  the  effect  of  replacing  [5]  by 


[B]  =  -3s  [I]  + 


resulting  in  the  modified  bound, 


—2(1  +  ar) 
(1  +  am) 


2(1  +  ar) 
-(1  +  am) 


A  f  < 


V  + 


Ax 


(78) 


(79) 


The  reason  we  opr  for  this  ad-hoc  bound  is  that  it  provides  a  simple  modification  of  the 
bound  of  Eq.  (77)  and  avoids  the  solution  of  a  quadratic  equation  at  each  time  step,  which 
would  be  necessitated  by  a  consistent  analysis.  As  demonstrated  in  the  next  subsection, 
this  bound  is  very  good  for  a  =  0  and  only  slightly  pessimistic  for  a  ^  0.  We  should  also 
point  out  that,  from  Eqs.  (Go)  and  (71)  and  from  the  fact  that  the  added  term  in  Eq.  (71) 
is  0(  At),  it  is  an  easy  matter  to  verify  that  the  ODQ  method  is  consistent  but  no  more 
than  first  order  accurate. 


D  Xuiwricnl  results:  Figs.  10  —  13  summarize  the  stability  and  accuracy  of  the  ODQ 
method  for  fiber  drawing.  For  a  discussion  of  the  nature  of  the  results  and  their  physical 
interpretation,  see  Ref.  [22].  The  results  displayed  in  these  figures  were  obtained  by  using 
the  bound  of  Eq.  (79)  to  compute  a  virtually  steady  solution,  which  was  saved  for  the 
stability  and  accuracy  studies.  The  “true”  stable  time  step,  labelled  A t  on  the  plots,  was 
obtained  by  interval  halving  from  a  bracket  obtained  by  incrementing  the  bound  of  Eq.  ( 79 ) 
by  a  sufficient  amount  to  obtain  an  unstable  step.  The  A t  reported  is  the  last  stable  time 
step  observed  at  the  point  where  the  bracket  had  been  reduced  to  a  small  enough  size  to 
make  no  difference  to  graphical  accuracy.  The  stable  time  step  thus  calculated  applies 
only  to  th<’  specific  steady  solution  involved,  but  the  bound  of  Eq.  (79)  has  been  found 
to  be  sufficient  for  stability  throughout  dynamic  processes  leading  to  steady  solutions. 
Tin'  convergence  plots  were  obtained  from  an  estimate  of  the  exact  solution  obtained  by 
Richardson  extrapolation  to  0fAx3).  Since  the  expected  accuracy  is  first  order,  and  the 
extrapolations  to  second  order  agree  very  closely  with  the  extrapolations  to  third  order, 
flu’  assumption  that  the  values  extrapolated  to  third  order  are  exact  with  respect  to  the 
raw.  first  order  results  seems  well-justified.  Note  that  the  plotted  errors  are  absolute  errors, 
and  the  relative  errors  are  quite  acceptable  with  ICO  grid  cells. 

There  are  two  major  points  to  be  made,  based  on  the  numerical  results:  First,  the 
stability  bound  of  Eq.  (79)  works  well,  and  second,  with  small  e.  acceptably  accurate 
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Fig.  10.  Mesh  Refinement  Studies  —  Corotntional  (a  =  0)  Case: 
(A)  Actual  At.  compared  At  predicted  by  the  stability  bound  of 
£(/.  f  70).  and  to  the  usual  wave-speed  bound,  in  the  for  the  £  =  1 
c; is**.  ( D )  Convergence  rate  of  maximum  values  of  key  quantities  to 
the  value  obtained  by  multi-mesh  extrapolation.  Slope  is  essentially 
-1.  indicating  a  first-order  method. 


results  arc  obtained  long  before  the  conventional  bound  of  At  <  Ax / v  becomes  useful 
on  very  tine  meshes;  thus,  the  bound  of  Eq.  (79)  is  necessary  to  obtain  solutions  in  the 
practical  operating  range  of  the  method.  Fig.  10  shows  the  corotational  case  with  a  large 
there  is  no  stiffness  in  this  case,  and  as  indicated  in  panel  (A),  the  conventional  bound 
becomes  adequate  on  meshes  of  SO  grid  points  or  more  and  could  be  used  generally  with 
a  modest  safety  factor.  The  expected  first  order  convergence  rate  with  grid  refinement  is 
confirmed  here  and  in  all  other  cases.  The  stability  situation  is  markedly  different  with 
small  £.  as  illustrated  in  Fig.  (11).  With  80  to  100  grid  points,  where  solution  accuracy 
is  quite  good.  th<'  stability  bound  of  Eq.  (79)  is  very  sharp  and  differs  substantially  from 
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Fix'.  IF  Mesh  Refinement  Studies  —  Corot  at  ional  (a  =  0)  Case: 
(A)  Actual  -At.  compared  A t  predicted  by  the  stability  bound  of 
E<[.  (79).  and  to  the  usual  wave-speed  bound,  m  the  for  thes  =  0.01 
rase.  (D  )  Convergence  rate  of  maximum  values  of  key  quantities  to 
the  value  obtained  by  multi-mesh  extrapolation.  Slope  is  essentially 
-1.  indicating  a  first -order  method. 


tin1  conventional  bound:  steady  solutions  could  not  be  obtained  on  such  grids  using  the 
conventional  bound. 

In  the  upper  convected  case  with  large  £  shown  in  Fig.  (12),  the  conventional  bound 
is  conservative,  and  the  bound  of  Eq.  (79)  is  slightly  more  so,  but  both  bounds  do  well. 
When  is  small  (Fig.  (13)),  the  conventional  bound  is  over  optimistic  until  the  grid  has 
G40  points.  The  bound  of  Eq.  (79)  is  slightly  conservative,  but  does  well  throughout  this 
range.  For  grids  finer  than  G40  points,  it  appears  that,  as  with  the  large  £  case,  the 
conventional  bound  is  becoming  a  viable  bound,  slightly  less  pessimistic  than  the  bound 
of  Eq.  (79):  however  for  very  fine  grids,  both  bounds  are  essentially  the  same.  From  a 
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Fig.  12.  Mesh  Refinement  Studies  —  Upper  Convected  (a  =  1 ) 

Case:  (A )  Actual  At.  compared  At  predicted  by  the  stability  bound 
of  Eq.  (70).  and  to  the  usual  wave-speed  bound,  in  the  for  the  £  =  1 
case.  (D)  Convergence  rate  of  maximum  values  of  key  quantities  to 
the  value  obtained  by  multi-mesh  extrapolation.  Slope  is  essentially 
-1.  indicating  a  first-order  method. 

practical  standpoint,  the  bound  of  Eq.  (79)  is  sufficient  throughout  the  range  of  realistic 
grids  and  is  always  relatively  sharp. 

5.  CONCLUSIONS 

Much  of  what  wo  know  about  the  spurt  phenomenon  in  Poiseuillc  flow  is  now  accessible 
to  analysis  but  was  first  uncovered  by  numerical  simulation.  We  hope  that  a  similar 
situation  will  come  about  in  the  analysis  of  step-strain  experiments.  Our  preliminary 
numerical  investigation  suggests  that  a  deeper  analysis  using  non-monotone  constitutive 
equations  would  be  worthwhile.  From  a  rheological  point  of  view,  it  is  important  to 
discover  whether  the  failure  to  satisfy  the  Lodge- Meissner  relation  is  an  indication  of 
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Fig.  13.  Mesh  Refinement  Studies  —  Upper  Convected  (a  =  1) 
Case:  (A)  Actual  At,  compared  At  predicted  by  the  stability 
bound  of  Eq.  (79),  and  to  the  usual  wave-speed  bound,  in  the 
for  the  £  =  0.01  case.  (B)  Convergence  rate  of  maximum  values  of 
key  quantities  to  the  value  obtained  by  multi-mesh  extrapolation. 
Slope  is  essentially  -1,  indicating  a  first-order  method. 


defect  in  the  constitutive  equation,  as  is  widely  believed,  or  the  failure  represents  a  success 
of  the  constitutive  equation  in  predicting  a  new  kind  of  material  response  to  step  strains 
that  can  occur  in  certain  (perhaps  extreme)  conditions.  Our  numerical  simulations  have 
made  the  connection  between  material  response  in  step  strain  and  the  spurt  phenomenon 
that  may  prove  important  in  resolving  this  rheological  controversy. 

From  a  technical  point  of  view,  the  major  conclusion  is  that  the  Newtonian  time 
scale  must  always  be  respected.  When  this  is  very  short  compared  to  the  dominant  relax¬ 
ation  time,  the  resulting  stiffness  affects  accuracy  alone  for  the  implicit-explicit  shear  flow 
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method.  This  stiffness  affects  stability  and  accuracy  for  the  ODQ  method.  In  either  case, 
A t  —  0(e)  is  often  called  for.  Fortunately,  however,  in  the  problems  we  have  solved  so 
far,  mesh  size  scales  of  0(e)  have  not  been  found  necessary;  in  fact,  the  major  technical 
problem  with  source  term  stiffness  is  that  quite  accurate  solutions  can  be  obtained  on  grids 
that  are  not  fine  enough  to  render  the  conventional  bound  on  At  useful. 

The  ODQ  method  seems  to  show  some  promise.  The  results  we  have  obtained  so  far 
lead  us  to  believe  that  further  study  of  the  ODQ  method  is  warranted. 
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